Draped Hypsometric DEM Workflow

# define the coordinates of Oroville Dam
oroville_df <- data.frame(
  Y = 39.534344,     # latitude
  X = -121.486103    # longitude
)

# convert to an sf point object
oroville_sf <- sf::st_as_sf(
  oroville_df,
  coords = c("X", "Y"),
  crs = 4326          # WGS84 geographic coordinate system
)

# load US state polygons
us_states <- maps::map_data("state")

# select only California
cali <- us_states[us_states$region == "california",]

# convert polygon coords to sf geometry
polygon <- sf::st_polygon(
  list(as.matrix(cali[, c("long", "lat")]))
)

# create an sf object from polygon
cali_sf <- sf::st_sf(
  id = "california",
  geometry = sf::st_sfc(polygon),
  crs = 4326
)

# download HUC8 polygons
huc8 <- nhdplusTools::get_huc(
  id = c("18020121", "18020122", "18020123"),
  type = "huc08"
) %>%
  st_transform(4326)   # ensure consistent CRS

# dissolve multiple HUC8s into one combined polygon
feather_dissolve <- huc8 %>%
  summarise(geometry = st_union(geometry)) %>%
  st_as_sf()

# compute bounding box
feather_bbox <- st_bbox(feather_dissolve)

# convert bbox to polygon
feather_bbox_sf <- st_as_sfc(feather_bbox)

# download NHD flowlines
nhd_flow <- nhdplusTools::get_nhdplus(
  AOI = feather_dissolve,
  realization = "flowline"
) %>%
  st_transform(4326)

# Add custom linewidths by stream order
# cartographic hierarchy
nhd_flow <- merge(
  nhd_flow,
  data.frame(
    streamorde = 1:6,
    width = rev(c(1, .8, .6, .4, .2, .1))
  )
)

# download lakes/ponds
nhd_wb <- nhdplusTools::get_waterbodies(
  AOI = feather_dissolve
) %>%
  st_transform(4326)

# download DEM tiles
elevation_data <- elevatr::get_elev_raster(
  locations = feather_dissolve,
  z = 10,                         # zoom level
  prj = st_crs(4326)$proj4string
)

# clip DEM to watershed
elevation_data <- raster::mask(
  elevation_data,
  feather_dissolve
)

# convert raster to terra object
r <- terra::rast(elevation_data)

# compute slope and aspect
slope  <- terra::terrain(r, "slope",  unit = "radians")
aspect <- terra::terrain(r, "aspect", unit = "radians")

# generate hillshade raster
hill <- hillshader::shade(
  slope,
  aspect,
  angle = 30,      # sun elevation
  direction = 270  # sun azimuth
)

names(hill) <- "shades"

# hcl.colors() — grayscale palette for hillshade
pal_greys <- hcl.colors(1000, "Grays")

# Rescale hillshade values to palette indices
index <- hill %>%
  mutate(
    index_col = scales::rescale(shades, to = c(1, length(pal_greys)))
  ) %>%
  mutate(index_col = round(index_col)) %>%
  pull(index_col)

# Map grayscale values to palette
vector_cols <- pal_greys[index]

# draw rasters
map_hypso <- ggplot() +
  tidyterra::geom_spatraster(
    data = hill,
    fill = vector_cols,
    maxcell = Inf,
    alpha = 1
  ) +
  tidyterra::geom_spatraster(
    data = r,
    maxcell = Inf,
    show.legend = FALSE
  ) +
  tidyterra::scale_fill_hypso_tint_c(
    limits = c(0, 2730),
    palette = "wiki-2.0_hypso",
    alpha = 0.6
  ) +
  theme_minimal()

# Add lakes/ponds
map_hypso_nhd <- map_hypso +
  ggplot2::geom_sf(
    data = nhd_wb[nhd_wb$ftype %in% "LakePond",],
    fill = "darkslategrey",
    color = NA
  ) +

  # Add flowlines with hierarchical widths
  ggplot2::geom_sf(
    data = nhd_flow,
    ggplot2::aes(linewidth = width),
    color = "darkslategrey", show.legend = F
  ) +
  scale_linewidth(range = c(0.05, .3))

# Build inset map
inset <- ggplot() +
  geom_sf(data = cali_sf, fill = "gray90") +
  geom_sf(data = feather_dissolve, fill = "grey20") +
  geom_sf(data = feather_bbox_sf, fill = NA, color = "black") +
  theme_void()

# compose final map
final_map <- cowplot::ggdraw() +
  cowplot::draw_plot(map_hypso_nhd) +
  cowplot::draw_plot(inset, x = .73, y = .65, width = 0.25, height = 0.25)